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Abstract 

Despite their importance in activated processes, transition-event durations — which are much 
shorter than first passage times — have not received a complete theoretical treatment. We therefore 
study the distribution pi,{t) of durations of transition events over a barrier in a one-dimensional 
system undergoing over-damped Langevin dynamics. We show that phU) is determined by a 
Fokker-Planck equation with absorbing boundary conditions, and obtain a number of results, 
including: (i) the analytic form of the asymptotic short-time behavior (t — > 0), which is universal 
and independent of the potential function; (ii) the first non-universal correction to the short-time 
behavior; (iii) following Gardiner a recursive formulation for calculating, exactly, all moments 
of Pi, based solely on the potential function — along with approximations for the distribution 
based on a small number of moments; and (iv) a high-barrier approximation to the long-time 
(t — > oo) behavior of pbit)- We also find that the mean event duration does not depend simply on 
the barrier-top frequency (curvature), but is sensitive to details of the potential. All of the analytic 
results are confirmed by transition-path-sampling simulations, implemented in a novel way. 
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I. INTRODUCTION 



Stochastic descriptions of dynamics have been invoked to describe protein folding 
protein dynamics j^, chemical isomerization [4], and chemical reactions 1^ among 
many other examples too numerous to list here. A Brownian particle moving in a one- 
dimensional bistablepotential typically provides a model for chemical and biological 
reaction systems |6tl3'|8|]. 

Basic aspects of the problem we address can be understood by examining Figs. [Hand 121 
where we show a trajectory for a Brownian particle moving in a one-dimensional double- 
well potential. There are two timescales of interest, see |9|]. One is the waiting time, or 
first passage time (FPT), which is the time the particle stays in one potential minimum 
before it goes to the other minimum. The other timescale is that for just climbing over 
the barrier separating the two minima, excluding the waiting time. We refer to this latter 
time as the "transition-event duration"; it has also been termed the "translocation time" 
in the context of membrane and pore traversal jlfl|. In the right hand graph of Fig.lll one 
transition-event extracted from the full time series on the left is shown at higher temporal 
resolution. See also Fig. |2| 

Theoretical analysis of the first passage time is largely a textbook subject now 111. [ill. 
In the simplest description, activated dynamics are modeled as Poisson processes, and 
the first passage time hence follow an exponential distribution. The mean FPT is given by 
Kramers' theory | ^Q], and depends on the height of the barrier as well as the curvatures 
at the minimum and at the barrier top. More complex models of the FPT distribution 
have also been studied HQ- 
Transition event durations have received focused attention only more recently 0, 
The event duration is important because it reflects the detailed dynamics 
of an activated transition. The statistics of event durations are of biological interest in 
many situations, such as transport in ion channels polymer translocation 



through a pore . Recently developed nanosecond and femtosecond-scale 

experimental techniques I20I. 2l|] have the potential to probe, directly, transition events. 

so fundamental to simulation techniques, such as path sampling 



This short timescale is a 
HE, 12112 



2^ 3I26 , 27 , 28 , 29 ]; see also (s^ 31 , 32 , 33 ,QSSl- We further note that 



many 



previous theoretical treatments based on optimization of the Onsager-Machlup action 
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discuss aspects of the transition event duration H-lsg , 



40 



41 



Also, the distribution of transition paths is directly described in the work of Dykman and 



colleagues ||4g 

The transition-event duration may be defined in a simple way In particular, Fig|2l 
shows that two positions are defined as the start and end points of the transition; usually 
they are on the different sides of a barrier. When a transition occurs, the event duration 
is the time interval between the last time the particle passes the start point and the first 
time it reaches the end point. Although arbitrary choices are inherently necessary for the 
start and end points, these do not appear to affect the basic physics. 

The probability distribution of transition-event durations was previously studied in a 
phenomenological way by Zuckerman and Woolf [4]. Hummer [51] gives an analytic 
formula for the mean transition-event duration for an arbitrary one-dimensional potential. 
Indirectly, Redner's study of the first passage time in an interval supplies important 
precedents for our work [52], as does Gardiner's book llj]- Two groups have recently 
discussed the time-reversal symmetry of the Oh distribution, albeit without attempting 
the detailed probe of the distribution itself lillla]/ which we pursue here. Other efforts 
directed at polymer translocation |117 , 11911 investigated a related but distinct problem, 
critically differing in boundary conditions; see Sec. II VI 

In this work, we first review the derivation of the probability distribution of the 
transition-event durations, pb{t), from the Fokker-Planck Equation (FPE) with particu- 
lar boundary conditions. We then obtain novel results. A recursive formula for all the 
moments of pb{t) is found, which permits accurate numerical approximations of pb for an 
arbitrary potential. The short-time behavior of pb is studied by path integral techniques, 
yielding universal behavior along with a potential-dependent correction. For a bistable 
potential with a high barrier (i.e. a "double-well"), the long time behavior of ph{t — > oo) 
will be discussed. 



II. BROWNIAN MOTION 

Our description of transition-events will be based on the traditional approach of a one 
dimensional "reaction coordinate" x coupled to a thermal bath Iz^. The analysis assumes 
over-damped Brownian dynamics, which we will variously address via a Langevin de- 
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scription, the associated Fokker-Planck equation and related path integral methods. In 
the following three subsections we introduce our notation, terminology, and the basic 
model. 

The direct object of our study, the distribution of transition-event durations, will be 
fully introduced in Section HVl 



A. Over-Damped Langevin Dynamics 

We consider Brownian, stochastic dynamics as governed by the over-damped Langevin 
equation for a generalized coordinate x in the presence of Gaussian white noise 1271] . 

^ = ^ + R(0, (1) 

at y 

where 

is the physical, systematic, conservative force acting on the particle, based on the potential 
energy U, and y is the friction constant. In this work, the noise, R{t), is taken to be Gaussian 
and white with zero mean and correlation 

{R{t)R{t'))={^y{t-t'). (3) 

Here fcg is Boltzmann's constant, and T is the temperature. 

Individual realizations of the noise in Eq.((lll generate stochastic trajectories x{t), which 
are routinely simulated numerically as described in Section Hill 



B. Fokker-Planck Equation 

To extract statistical information on trajectories one generally turns from the Langevin 
equation to the associated Fokker-Planck equation (FPE). The relation is discussed in 
standard monographs; see, for example [6]. The FPE describes the average behavior of a 
statistical ensemble of trajectories x{t). 

In one dimension the time evolution of the probability density function P{x, t), for the 
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coordinate x is assumed to be described by the Fokker-Planck equation 



dP{x,t) 

dt 



P(x, , 



(4) 



In Eq.(|4|l, D^^\x) > is the diffusion coefficient and D^^\x), is the drift coefficient. When 
D*^' = k^T ly = D, with D the diffusion constant, and D^^\x) = F{x)/y, this Fokker-Planck 
equation embodies the over-damped Langevin dynamics of Eq.((lll with noise correlation 
satisfying Eq.(|3l( |^. In this simplified case we obtain 



dP{x, t) _ d i 
dt dx 1 



dU*(x) d 

— H 

dx dx 



F{x,t)\ , 



(5) 



where ll*{x) = U{x)/kBT is the dimensionless physical potential. The solution of the 
Fokker-Planck equation under suitable initial and boundary conditions will allow statis- 
tical information to be extracted. 

Eq.© clearly expresses conservation of probability. The total current associated with 
the stochastic variable x is given by |l| 



}{x,t) = -D 



dW{x) d 

— H 

dx dx 



(6) 



where the first term represents the drift (systematic) contribution, while the second term 
is the diffusion contribution. 



C. Path Integral Approach 



Path integral methods provide a useful tool and a different perspective for the study 

, |5J, 55]. Trajectories x(f) are directly considered, rather than 



of Brownian motion 



53 



their average as in the EPF. The relative probability of a trajectory, connects the 

path integral back to the original Langevin equation ((D). It can be shown that the relative 



probability for a Brownian particle described 
specific path x(t), with x the time, is given by 



jy the Langevin equation (HJ to follow a 



W[x(-)] = exp 



L[x(T)]dT 



(7) 
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where the effective Lagrangian can be expressed as 

, (dx IDdF 

Then the optimal (classical) path Xc{t) can be found in the standard fashion from the 
stationarity of the exponent "action". 



6 r L[x(T)]rfT = 0, (9) 

with x(io) = a and x{t) = x. Integrating over the probability densities of individual paths, 
the propagator 

G{x,t\a,to)= r W[x(T)]d[x(T)] , (10) 

determines the probability of the particle arriving at position x = x{t) at time t given that 
it started at a = x(io). With a suitable definition of the measure in the integration in (ITO|) 
the propagator is equivalent to the "principal solution" or Green function solution of the 
associated Fokker-Planck equation |53l.l55|1 . 

When the diffusion coefficient D is small, the major contribution to the propagator will 
come from the paths very close to the optimal path Xc{t). So, for D — > 0, the simplest 
approximation retains only the contribution from the optimal path L53.1 , and is of the 
form 

G(x, to) = m exp L[Xc{T)]dT^ , (11) 

where K{t) is the normalizing factor to ensure the propagator satisfies 



G(x, t\a, to)dx = 1 . (12) 



To include small fluctuations around the classical path, one typically invokes a 
quadratic approximation, in which deviations to second order are retained in the effective 
Lagrangian. Writing 

x(t) = x,(t) + 6(t), (13) 
with 6{to) = 6{t) = 0, and neglecting terms in the Lagrangian higher than the second order 
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in 6, one finds 



G{x,t\a,to) = exp|^^^ Fdx'^(^4nD[2E + 4DV{a)]H2E + 4DV{x)]^ {2E + 4DV{x'))-'^dx' 



X exp 



Here the effective potential, V{x), (distinct from the physical potential U) is given by 

, F(x)2 1 dF 
V(x) = + ■ (15) 

The constant E in Eq. ((T4)l is the "energy" of the particle, appearing as an integration con- 
stant in the first integral of the effective equation of motion following from the Lagrangian 
in Eq.© and the extremization in Eq.©. The first integral yields the equation of motion, 

dx 1 

^ = {2E + 4Dy[x,(T)]}5, (16) 

which is the equation for a classical particle of mass | in the presence of a conservative 
potential -2Dy(x) (sS]. Larger values of E give optimal solutions for increasingly more 
rapid events. 

The quadratic approximation is justified if the path probability decreases sufficiently 
rapidly with increasing variation of the path from the optimal one | Ssl. This is expected 
in the limit of weak diffusion (D — > 0) in analogy with the semi-classical approximation 
in quantum mechanics with ^ — > 0; see,e.g., [53]. 

The path integral approach provides some interesting insights for the features of Brow- 
nian dynamics of concern in this paper. We will return to this description below. 

III. SIMULATION 

All of our key analytic results to be discussed below have been confirmed via numerical 
simulation. Here our simulation approaches are briefly described. 



(14) 
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A. Brute Force Simulation 



Standard simulations of the over-damped Langevin Eq.Q employ a simple first-order 
scheme with fixed time step M, such that 



Xj = xijAt), i = 0,1,2... 



(17) 




F{xj) 



At + Ax 



(18) 



y 



R ■ 



Consistent with Eq.(|3ll the thermal fluctuation (noise increment) Axr is chosen from a 
Gaussian distribution of zero mean and variance 



However, as is well-known, this direct approach proves inadequate to simulate rare 
events, even in one dimension. A program running on a single CPU can provide an 
ensemble of transition trajectories (with thousands of transition-events) only for low 
barrier height. For high barriers, the waiting time between successful events will become 
unacceptably long. Therefore, we employed a path-sampling method for simulations, 
which we now describe. 

B. Path Sampling 

The problem with direct "brute force" simulations is that the waiting time between 

n 

events grows exponentially with barrier height II ID . Our interest here, moreover, is to 
obtain a statistically well-sampled ensemble of transitions. In practice, for any given 
model, we require thousands of events. 

To generate a sufficient quantity of transition-events, we turn to a Monte Carlo path- 
sampling approach. The approach has its roots in Path Integral Monte Carlo for quantum 
systems Isdls^ . but Pratt provided an important advance in recognizing the analogous 
application in classical and, particularly, chemical systems [22]. 

Pratt's approach has recently been taken up with some vigor by Chandler and co- 




(19) 
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workers Related work was presented by Zimmer and Paniconi I59L l60|1 . An 

independent path-sampling approach was developed by Zuckerman and Woolf Iz^ . 
building on work by Ottinger 16 ill . 

The basic idea of path-sampling is simple: focus computer time on the rare transition 
events of interest. Fig. fright), rather than on the waiting time between events, which can 
be longer by many orders of magnitude, Fig.lHleft). In a statistical mechanics context, the 
probability of apath can be computed. Hence "trial" paths can be included by means of 
re-weighting |8|] or by a Metropolis criterion i24\ . 

Here we primarily follow Pratt's approach to path sampling, which is based on two 
facts: (i) Path (i.e., trajectory) probabilities are readily computed for stochastic processes, 
so that trajectories may be viewed as Nxd dimensional equilibrium "objects", when there 
are N time steps and d spatial dimensions, (ii) Wherever equilibrium probabilities can be 
computed for all such "objects" in a space. Metropolis sampling can be performed. 

As in any Metropolis simulation, we require that detailed balance is satisfied. That is, 
for arbitrary paths i and j with equilibrium probabilities P and overall transition rates T, 
we require 

PpathUni ^ ;■) = PpathUm ^ . (20) 

The rate T is decomposed into the usual product of the generating (gen) and acceptance 



(ace) components 11241. 16211 . which are proportional to the conditional probability for gener- 
ating and accepting the trial path starting from i. Then trial moves should be accepted 
with probability min[l, J^], where 

^ ^ accji i) ^ Ppathii) X genii i) 
accij i) Ppathii) X genii j) ' 

All paths in our ensemble will have the same total number of steps N, so that the 
probability of two paths can be compared via Eq.(|2H). We will typically choose N to be 
much bigger than {t)b/At, with {t)h the mean time for transition events, so that the "full 
shape" of the distribution is sampled. Thus, transition events will typically constitute 
only part of N-step trajectories in our sample of paths. On the other hand, intentionally 
selecting N smaller than {t)i,/At when necessary allows us to focus on the short-time 
behavior of pi,. 
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To proceed, we must establish the equilibrium and generating probabilities in Eq. (|2Hl . 
The "equilibrium" probability Ppath of the N-step path from a to is the product of the 
equilibrium probability for the initial point and all subsequent single-step transition prob- 
abilities consistent with Eqs.dTSb and (IT9b . We further restrict our ensemble to "successful" 
paths containing transition events by formally introducing a projection operator 6. Thus 
we have 

Ppathiia, xi, Xn}) oc exp[-U*{a)] x [U^^-'^p{xi, x;+i; U*)] x 0({x,}) . (22) 

The single-step transition probability corresponding to Eqs. (IT8|) and (IT9|) is a Gaussian 
density, namely 



p{Xi, x/+i; LT*) = — exp 
^/2no 



1 [ [x,,,-x,-l(ff)(2DA0f' 



2(72 



(23) 



where . ^ 



. If the particle returns to the left boundary a before arriving at the 
right boundary b, 6 = 0; otherwise 6 = 1. 

We employ a novel path generating procedure designed to focus simulation effort 
on the distribution pb{t) of interest and, as necessary, on the rarest trial events. Our 
path generation strategy is closely related to a non-Metropolis re-weighting procedure 
previously considered by Zuckerman and Woolf [8]. Specifically, a trial path is built 
up "from scratch", but based on the average behavior of the previous path. From the 
previous path, which starts from a and arrives at b after Mom {Mom < N) steps without 
being absorbed at a, we can calculate the average velocity over total time Moid{2DAt) as 

This will be the "target speed" of the new trial path. This is extremely useful when 
studying the fastest events, whose transition-event durations are much shorter than {t)h. 
To generate a new path, we linearly bias the particle from atob using 

xy+i = Xj + (viMoidWDAt) + Axr , (25) 
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where Axr has been defined in Section PlI Al following EqdTst. Eq. JZSb may be compared 
to the unbiased form dlSt . The linear bias in JZSt is motivated by the quasi-ballistic quality 
of the fastest transition-events deriving from Eq.(IT6b in the limit E » DVmaxr where Vmax is 
the maximum of V{x) defined by Eq.dTst. Thus, on the new path, the particle moves with 
a constant drift (bias) velocity, as if the force were constant, and ordinary noise. Note that 
for the new trajectory generated by Eq.dZSb, the new value, M„ei„, can be larger or smaller 
than Mold- Once it arrives at h, we remove the bias and allow the particle to move for the 
remainder of N steps as governed by unbiased Brownian motion Eq.dTst. As noted above, 
all the paths must contain the same number of steps for probabilities to be compared in 
our Metropolis procedure using Eq. (|2T|l 

The generating probability {gen) for our procedure is the conditional probability with 
which we choose the new path, given the old one (with its average speed), namely, 

M„™-1 N-1 

gen{old ^ new) = ]^ g(x„ x,+i; u(Mom)) ]^ p{Xj,Xj+i;U*) , (26) 

!=0 j=M„nv 

where 

g{xi, Xi+i, v{Moid)) = exp I — 1 . (27) 

This generating method is tailored to the potentials and boundary conditions we study 
in this paper, so that the Brownian particle will not be trapped in any position between 
the two absorbing walls. 

By substituting Eqs. (|22)l ,(l23b,(l26t and (|27|) into Eq.(|2T|) we arrive at the acceptance 
criterion for our generating procedure, namely 

R — r-f — — — , (28) 

given an old transition path (fl, 1/1,1/2 Vmm'Vmm+i Vn) and a trial transition path 

{a, Xi, X2 ^M,„,„/ ^M„„„+i ^n)- 

To ensure the correct behavior of our procedure and code, path-sampling results were 
carefully checked against direct simulation, using Eq.JTSt, in a number of cases. In the 
following sections of the paper, all the simulation results employ the path-sampling 
method just described. 
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We also checked that our path-sampling simulations greatly exceeded the correlation 
time resulting from our use of the "old" average velocity in Eq. (|24)l . This resulted in 
negligible statistical uncertainty, as can be gauged from the smoothness of the data in all 
path-sampling figures. 



IV. TRANSITION-EVENTS 



A. Distribution of Transition-Event Duration 



The distribution of transition-event durations, pb{t), for a Brownian particle confined 
to one dimension can be found by solving the Fokker-Planck equation using suitable 
boundary conditions, as we now describe lll^ . 

During the entire transition process the particle must move between a and b, which 
means that only trajectories that stay completely within the interval are considered, i.e. 
a < x{t) < b during the entire event. To eliminate the extraneous trajectories, absorbing 
walls must be put at the start and end points, ll^ . so that 

P{a,t) = 

P{b,t) = 0. (29) 

As recently stressed by Berezhkovskii |lo|, et al., the dual absorbing boundary conditions 
distinguish the event duration as a "conditional first passage time", rather than the usual 
unconditional time associated with the Kramers' pro blem I9L IIOII . This contrasts with 
several previous studies of polymer translocation 

One releases particles very close to the left absorbing wall at i = 0, so that the initial 
condition is 

P(x,0) = 6[x-(fl + e)], (30) 

with e — > 0-I-. Then the current at the right absorbing wall will determine the distribution, 
Pb, of durations according to 

pf,(0 cc lim J{b,t), (31) 
with the currents given in Eq.(|6ll. Following Gardiner's work 111], let 7if,(fl -I- e\t) equal the 
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probability that a particle, released ata + e, is absorbed at the right absorbing wall during 
< T < i. It is easy to see that 



r m- 

Jo 



Tibia + e\t) = I }{b, T)dT . (32) 
If we define 

}{b, T)dT , (33) 



this "splitting probability" can be used to normalize p\,{^) in Eq.dSTl) according to 

pdt) = lim = lim — — . (34) 

^-0+ Jo Jib,^)dT ^-o+nfc(fl + e) 

We note that the splitting probabilities are time-independent and follow directly from the 
potential U* according to jlilsil: 

f''exp[U*{x')]dx' 

n«W = ^ 

I exp[Lr*(x')]^^^' 

rexp[Lr*(x')]^^^' 
Ubix) = = 1 - a(x) . (35) 

JT exp[U*{x')]dx' 

Hence, to find the distribution of the transition-event durations, pbit), one must solve 
the Fokker-Planck equation (jSj with the initial condition (l30t and absorbing boundary 
conditions (l29t . The current, }{b, t), can be found from Eq.©, which can then be combined 
with the splitting probabilities to find the normalized distribution of transition-event 
durations, pbit). 

B. Examples: Free Diffusion and Linear Potential 

The solution of the Fokker-Planck equation can be formally expressed, in standard fash- 
ion, in terms of the eigenvalues and eigenf unctions of a time independent equation 
The solution can be written in the form 

P{x,t) = Y^A„pn{x)e-''"', (36) 
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where the eigenvalues {A„} are non-negative and, based on Eq.©, the eigenfunctions 
satisfy 

( d \dW{x)} . <f\ , , , , , ^^^^ 



dx 



Eq.(l37b with boundary conditions (l29t determine the functions {pn(^)}/ while the constants 
{A,,} are found from the initial condition (l30t . The distribution pb(i) follows from Eqs. (l34|l 
and (|35)l . As examples, we determine p;,(i) for a few special potentials Lr*(x). This will 
reveal some interesting features. We note that the linear potential, of which free diffusion 
is a special case, previously was studied by Lubensky and Nelson [14], although without 
numerical simulations. 

Free Diffusion. Even in the absence of a true barrier, the event duration is still well 
defined by the formalism above, and this simple case acts as a useful reference. We 
therefore first consider free diffusion, with U*{x) = and a = 0,b = L. 

The solution of Eq. ()37)l can easily be found, and the result can be formally expressed as 



P^it) = 2D 



oo / ■ \2 

;=1 \ / 



exp 



"IT 



Dt\ 



(38) 



Notice that at long times the decay is exponential and dominated by the lowest eigenvalue. 

The right-hand side of Eq. (|38)l is well behaved for long time, but is not useful for 
^ — > 0. We can re-cast the result in a format useful at short times by using the Poisson sum 
formula 

oo oo ^oo 

X/(")=Z I fix)exp{-2inix)dx, (39) 



for function /. We then find an alternative representation 



2L 



y 



(2; + lyi^ - IDt 
Wt 





r {2i + lfL^^ 


exp 


ADt 



(40) 



which can be used to extract the behavior as i — > 0, namely pf\t) ~ t exp[-L^/(4Di)] 



We note that Eq.(l40t can also be derived using an image method, as described below in 
Section HVEH 

Linear Potential. The solution for the linear potential, U*{x) = kx, which corresponds 
to a constant drift velocity, can also be formally written in terms of an eigenfunction 
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expansion 



sinh(fcL/2) 
kL/2 



exp 



[{kL/2f + fn^] 
L2 



Dt 



(41) 



Comparing with Eq. (|38|l , the result can be written 



smh{kL/2) 
Ej2 



exp 



L2 



Dt 



(42) 



In the left graph of Fig. |3l we show path-sampling simulation results following Sec- 
tion UnBl for pfc for free diffusion with U* = 0, L = 1.0 and U* = 0, L = 2.0. They are 
compared with the numerical evaluations of Eq.JSSt. The path-sampling simulations and 
numerical results from the eigenfunction expansions match very well. We changed the 
units of the vertical and horizontal axes, so that all the curves of pf^ will not depend on 
the width L, and the generic behavior is highlighted. 

In the right graph of Fig. |3lwe show path-sampling simulation results for the a series 
of parameterizations of linear potential: U* = 4.0x, L = 2.0, U* = 8.0x, L = 1.0, U* = 9.0x, 
L = 2.0, and U* = 18. Ox, L = 1.0, which are compared to numerical evaluation of Eq. ()4T|l . 
The simulation and numerical results again match very well. We again scaled the axes 
to emphasize that the shape of ph only depends on the value of kL, which is essentially 
the potential energy difference between the start and end points. These exercises add 
confidence to the path-sampling methods used here. 



C. Approximate Solution for Inverted Parabolic Potential 

As a first investigation of a more realistic potential, we employ a crude representation 
of absorbing boundary conditions. In Fig. HI inverted parabolic potentials are shown, 
one with open boundary conditions and the other with two absorbing walls(Lr* — > -oo; 
see, e.g., |ll|]). When the "barrier" is high i.e., when Lr*(0) » U*{a = -1) = U*{b = 1), 
a particle exiting the region a < x < b, has a small likelihood of returning with open 
boundary conditions, because of the rapidly increasing "downhill" forces external to 
the region. Thus, as long as there is a sufficiently high barrier, one might conclude the 
solution for open boundary conditions will be a good approximation for an inverted 
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parabolic potential with two absorbing walls. We now investigate this approximation. 

With open boundary conditions, the exact solution of the Fokker-Planck equation for 
an inverted parabolic potential, U* = -^ax^, is well known 

a ( a\xexp(-aDt) - af-] 

2.[1 -exp(-2«D0] -P|- 2n - exp(-2«D0] |-P<-"°')- <«) 

which satisfies the initial condition 



P{x,0) = 6{x-a). (44) 

For the same potential with absorbing walls at - W and W, we approximate the current 
from Eqs.(l43t and ^ with a = - W and b = W, 

I 4 y n) sinh(aDt/2) JsmhiaDt) 



For normalization we will need 



N= f ]{Ht')dt' , (46) 
Jo 

which is the total probability passing to the right of x = W. Under the influence of this 
inverted parabolic potential, this probability will not pile up but will flow toward x — > oo. 

f exp{-}:ax^)dx erfcfWVoT^) 
N = n.(-W) = I = ^^t^^^ V^/^^ , (47) 

from which one obtains the approximation 

p,(0 . I 5^1^ 1 ( M exp[-(^WV2)coth(.D^/2)] 

U[l-erf(WV«72)]J \ V 7i/ sinh(aDi/2) Vsinh(«DO 

where erf(x) = 1 - erfc(x) = exp{-z^)dz Is^ . In Fig.|5lwe compare the results from 

direct simulation and from Eq. (l48|) for inverted parabolic potentials with different heights. 
In the simulations the two absorbing walls are placed at a = -1.0 and b = 1.0; then the 
height of the barrier is given hy a/2. 
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As expected, this approximation improves with increasing barrier height. 



D. Moments of the Distribution of Transition-Event Durations 



1. Recursive Formula 

In studying a distribution, it is natural to investigate its moments. Gardiner 111 
provides an expression for the first moment of pb{t) (i.e., the mean time), which is also 
given by Hummer Isil] . Here we derive a recursive formula for all moments of ph{t), as 
suggested by Gardiner ll ]. 

Following Gardiner, we define gh{x, t) as the total probability that the particle is ab- 
sorbed at b after time t, given that it is released at position x at ^ = 0. Thus 



J{b, T)dT, 



(49) 



and we have the initial condition 



P(x',0) = 6(x' -x), 



(50) 



The limiting cases for gi,{x, t) are 



gb{x,t = 0) = Ubix) 
gh{x,t = oo) = 0, 



(51) 



where Tlh{x) is defined in Eq.(l35t . 

The nth moment, T„, of the exit time distribution for particles released at arbitrary 
a <x <h can be calculated from gh{x, t) according to 

T^{h,x) = - r r| '^^ = T^^^ r i"~'Sb{x,t)dt, (52) 

Jo dt[gh{x,0)\ Tlbix) Jo 



so that T„{h,a) is the nth moment of pf,. Gardiner (1|] shows that gh{x,t) satisfies the 
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backward Fokker-Planck equation 



dx / dx 



dx^ 



dghjx, t) 
dt 



(53) 



Multiplying by ni" on both sides and integrating with respect to t yields 



-D 



dU* 
dx 



Jo 



^n-l dgb{x,t) 

dx 



dt 



+ D 



J 

Jo 



n i 



dt 



Jo 



n-l dgh{x,t) 
dt 



dt. (54) 



Now the right side can be integrated by parts to find 

,n-ldgb{x,t) 



n I 

Jo 



dt 



dt = -nIlb{x)Tn-iib,x) , 



and with Eq. 



-D 



dU*> 


j dy{x) 


dx j 


dx 



+ D- 



dx^ 



-nnb{x)Tn-iib, x), 



where 



y{x) = nh{x)Tn{b,x) . 
The boundary conditions on y{x) are 



(55) 



(56) 



(57) 



y{b) = y{a) = 0. 



(58) 



One way to solve equations like (l56b uses Green's functions Ig^ . The function that 
satisfies the homogeneous equation corresponding to Eq. (l56b with boundary condition 
(|58)l at a is, in fact, Tla{x); correspondingly at b, it is Ylb{x). Using these solutions, one 
obtains a recursive formula for all the moments 



Tr,{b,x) = ^^P[^*(^')]^^'}{n^ J exp[-U*{x')]nl{x')T„.,{b,x')dx' 



+ 



■L7*(x')]n„(x')n,(x')r„_i(fc,x')rfx' 



(59) 
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Our main interest is in the moments of pb, namely, 

Tn{b,a) = -^l^ exp[Lr*(x)]rfx| 1^ exp[-Lr*(x)]n„(x)nfo(x)r„_i(&,x)dx| . (60) 



Given the moments according to Eq.(l60t, the distribution of transition-event durations, 
pb{t), can be reconstructed numerically, at least for a fixed range of t. 



2. Lowest Eigenvalue 



From Eqs.© and (l36t one knows that pb(f) can be written in the series 



pb 



(61) 



where the eigenvalues defined by Eq. (|37|l satisfy < Ai < A2 < A3 < .... The eigenvalues, 
in particular Ai, can be found via direct numerical solution of Eq. (|37)l . Here we show 
an alternative based on integrations involving the potential. The first eigenvalue Ai can 
be expressed in terms of the high-order moments because of asymptotically exponential 
behavior. When n » 1, 



T„{b,a) 



-f 

Jo 



t"pb{t)dt 



^ A" 



1 + , Cs/Ai^""^ 



Ci VA 



^ Ci Ia3 



-I- ... 



(62) 



where T{n) is the Gamma function. The lowest eigenvalue can then be estimated from a 
ratio of high moments; for example. 



Ai - lim 



T„ib,a) 
Tn+i{b,a) 



(n + l) 



(63) 



and from Eq. (l62|) , the constant can be determined according to 



n+l 



Ci = lim 



Tn{b,a)A'l 
T{n + 1) 



(64) 
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Recalling that the moments can be constructed via successive integration, Eq.(l63b pro- 
vides a way to estimate the first eigenvalue in Eq.dSTI). In Section FlV Fl Eq. (l63t will be 
used together with simulations to check an approximate analytic result for the leading 
eigenvalue in a representative case. 



3. First Moment 



For n = 1, using To{b,a) = 1, Eq. (|60|l yields the first moment of the distribution of 
transition-event durations, 

Ti{b,a) = -^1^ exp[L[*(x)]dx| I exp[-LJ*(x)]n«(x)n;,(x)rfx| . (65) 



We can immediately evaluate Ti for the simple potentials. For free diffusion with U* = 
and a = 0,b = L, 

2DT? 1 

(66) 



2DTJ 1 
IT " 3 



With U* =kx,a = 0,b = L, 



2DTi 



L2 (kiy 



itLcoth(yj-2 



(67) 



For an inverted parabolic potential U* = H(l - ■^), where the curvature a = and 
a = -W, b = W, we can find an approximation of Ti. When H » 1, by using the method 
of steepest descents. 



and 



exp[Lr' 

w 



,2n 

{x)]dx ~ exp(H) -I / — , 



a 



n„(x)n;,(x) 



(68) 



(69) 



Then 



2DTi 



1^ r I 



1 - 



erf 



VHx^ 
W J 



dx 



■2 /-Vh 



exp(y2){l - [erf(y)]2}dy 
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H 



exp(y2){l - [erf(y)]2}dy + exY>{y^){l - [eri{y)f]dy 



^ ^[1.27 + log(H)] = -[1.27 + log(H)] . (70) 

Eq. (|70)l provides a good approximation for Ti for the inverted parabolic potential with 
high barrier. It also gives a rough estimate of Ti for a "single-bump" barrier with height 
H and width 2W. 

Naively, one might guess Ti should simply be proportional to the effective frequency 
for the barrier top, namely, the inverse curvature a~^. However, this intuition falls short 
in two respects. First, the logarithmic term in Eq.dTOt is dominant for large barriers, 
even for this simplest purely parabolic potential. Further, in extensive numerical work, 
for a double-well potential, we have seen unambiguously that the mean event duration 
is sensitive to details of the potential far from the barrier top (data not shown). This 
sensitivity can be traced to the dependence of the optimal "speed", Eq.(IT6t, on details of 
the potential. To give an extreme example, if there were a second barrier and minimum 
in the potential, then Ti would have to include the Kramers' time for the second barrier. 

4. Reconstruction of pi, from Moments 

Reconstructing a function approximately from a finite number of moments has been 
studied, e.g., by maximum entropy method 1671. 16811 . continued fraction approach 
Tol . and Talenti method ItiIIt^I and perhaps other techniques. Here we follow Hon 
and Wei's work 1721 to reconstruct the density pb{t) in a similar way. 

First one builds up an orthonormal set of basis functions ^j{t), which are polynomials. 



xl>j{t) = Y,Cjnt\ (71) 



n=0 



by using the standard Gram-Schmidt orthonormalization technique The polynomi- 
als satisfy 

I iPjit)4>jit)wit)dt = 6jj^, (72) 
Jo 
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with respect to a weight function zv{t), which is tailored to our problem with the choice 
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wit) = exp(-AiO exp j (1 + rS) , (73) 

where L = b - a. For this weight factor, when t ^ oo, iv{t) ~ exp(-AiO/ and when ^ — > 0, 
iv{t) ~ exp^-^^i"5. These forms represent the long-time and short- time behaviors for 
Pb{t) as we will show in Section llV Fl and ()IV E 4|l below. Thus we built in all the information 
we know about pb{t) in this weight factor. Notice that Ai can be found by following the 
scheme in Section llV D 21 or via other numerical methods. 

Following the usual Gram-Schmidt procedure, one builds up 

m 

pUt) = w{t)Y^aj^Pj{t), (74) 

;=0 

which will be an approximation for ph{t) of "order" m. If we incorporate moments of pi, 
by setting coefficients according to 



do = 



J-^oo poo 
g :po{t)w{t)dt Coo Jo Mt)dt 

f tiipj{t)w{t)dt 



«; = r./ i (75) 



then pni{t) will reproduce the first m moments T,, (0 < n < m), i.e., 

T„ = I t"pi,{t)dt ~ I ep„Xt)dt . (76) 
Jo Jo 

By using the first five moments, we reconstruct the distribution of transition-event 
durations for several different double-well potentials. Two of the results are shown in 
Fig. |6l They match well with the simulations except for the long time tail as seen in the 
semi-log plot. However, the event probability in that region is quite small. 
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E. Short Time Behavior 



Beyond the moments, it is of interest to study the asymptotic behavior of the event- 
duration distribution, in both short and long time limits. We first analyze the i — > 
behavior, using exact methods. 



1. Short Time Behavior for the Green Function with Open Boundary Conditions 



From the perspective of the path integral, introduced in Section HTCl if {t - to) is short, 
the velocity on the optimal trajectory ^ will be large. This implies a large "energy" E in 
Eq. ((T4|l . We therefore assume 2E » 4DV and obtain the corresponding Green function 
for short time t with open boundary conditions. 



G{x,t\a,to) ~ exp 



Fdx' 



^hT J a ' "" I ^4nD{t - to) 
X exp[-V{t-to)], 



exp 



(x - a)^ 



4D(i - to) 



(77) 



where V is the average effective potential between a and x, 



V 



1 r 

~ x-aX 



V{x')dx' 



(78) 



We have used the quadratic approximation, which is expected to be reasonable if the 
diffusion coefficient D is small 153(1 . Also, because we are further restricting our analysis 
to short-time behavior, the important paths will be close to the optimal one, which should 
improve the approximation. 

If we define Go as the Green function for free diffusion with open boundary conditions 

fl, 

{x - aY 



Go(x, t\a, to) = 



^A:nD{t - to) 



exp 



4D(i - to) 



(79) 



Eq. ((77|l can be expressed as 



G{x, t\a, to) ~ Go(x, t\a, to) exp ^ f j exp[-y(^ - ^o)] , 



(80) 



which will be useful when we discuss the early time behavior with alternative boundary 
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conditions below in Section PVE 21 



2. Short Time Behavior for the Current with Absorbing Boundary Conditions 

In this subsection we need to retrace the path integral method in order to study 
transition-event durations, as required, with two absorbing walls. 

As shown in Fig. |7l we wish to calculate the path integral from the start point '+' at 
X = 0+ to the end point 'o' at x = L- during the time interval {t - to). There are two 
absorbing walls at position x = and x = L, with some arbitrary potential between them. 
If the position of the start point is ei and the end point is x, then the Green function G"'"^ 
with absorbing boundary conditions can be expressed as a sum over Green functions G 
for open boundary conditions Is^l 



C^'ix, t\ei, to) = G(2;L + x, t\ei, to) - ^ G(2;L - x, t\e,, to) . 



(81) 



The construction is shown schematically in Fig.|7l 

As shown above, to determine the distribution of event durations, we first calculate 
the current at x, then take the limits ei — > and x — > L, 



/*(L,f|0,io) = -2D^ 



dG{x, t\0, to) 



dx 



x=(2/+l)L 



(82) 



From Eq. (|80|l , when {t - to) is small 



dG{x,t\0,to) 



dx 



dGoix,t\0+,to) 



x=(2i+l)L 



dx 



i-=(2/+l)L 



X exp 



1 r 



.(2;+l)L 



Fdx' 



xexp[-Vj{t-to)], 



(83) 



where the symmetry of the periodically continued potential has been used; see Fig.lZl 
The periodicity implies further simplifications, including 



Fdx' = 
Jo Jo 



Fdx' 
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(2; + l)L Jo 



.(2;+l)L 



V{x')dx' = - V{x')dx' = Vo . 
^ Jo 



If we define the free diffusion current with absorbing boundary conditions as 



/f(L,i|0,io) = -2D^ 



(9Go(x, i|0, 



dx 



.T=(2/+1)L 



(84) 



(85) 



Eq.(l82b can be written 



f%L, m, to) « JtiL, t\0, to) exp ^ Fdx'\ exp[-Vo{t - k)] , 



(86) 



or 



where 



f'{b,t\a,to)^JT{b,thk)exp 



abs/ 



,2/crT 



Fdx'\exp[-V{t-to)] 



V=-^ f V{x')dx' . 
We can now estimate the short time behavior of the normalized current. 



(87) 



(88) 



3. Short Time behavior for the distribution of the transition-event durations 



Combining Eq. (|34|l and (|87|l , we find that the short-time behavior for the distribution 
of transition-event durations is given by 



p„(i^O)^pf(Oexp 



' 1 



Fdx'\exp{-Vt) 



lim 



nlia + e) 



(89) 



where n||(x) is the splitting probability for free diffusion and is the distribution for 
free-diffusion with absorbing boundary conditions. From Eq. (l35|) 



n°(fl + e) exp[U*{x')]dx' 
e^Q+ Tibia + e) exp[U*{a)]{b - a) 



(90) 
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Combining Eq.^ and and using F{x) = hT[-^), the normalized distribution 
becomes 



pb{t ^ 0) 



f\xp[U*{x')]dx' 
exp[U*{a)](b - a) 



exp 



U*{b) - Wia) 



exp(-W), 



(91) 



revealing corrections to the free diffusion result due to the potential. 

In Fig. Is] we compare the results from a path-sampling simulation as described in 
Sec. nil Bl and our final result Eq.(l9Tt for two double-well potentials of varying barrier 
height. The simulations and the analytic results of Eq.(l9Tt are in good agreement at 
sufficiently early times, although at the earliest times the simulations reveal degrading 
statistics. 



4. Comment on the Short Time Behavior of (t) 

For the free diffusion problem with absorbing walls, the path integral method in 
Section liV E 21 can give the solution. Redner [52] shows that an image method, similar to 
that used in electromagnetism, will give the same result 

If Dt «; L^, 

p°(0~riexp(-^J. (93) 

Eq. (|93)l shows the short time behavior of p\, for free diffusion with absorbing boundary 
conditions. Combining it with Eq. (|9H) , one can find the short time behavior in the presence 
of the potential. 



(2; + lyu - IDt 
Wt 



exp 



(2; + l)^ L 

4Di 



2r2l 



(92) 



F. Long Time Behavior for Double- Well Potential 

From Eqs.© and ()36)l we know that the long time behavior of ph{t) will be determined 
primarily by the first eigenvalue Ai. Here we use a perturbative approach to obtain 
an approximation for Ai for a double-well potential with high barrier. We confirm, by 
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direct numerical calculation, the validity of using perturbation theory. We also perform 
path-sampling simulations to check the accuracy of our final approximation. 

By a variable transformation, the one dimensional Fokker-Planck equation can be 
transformed to a Schrodinger-like equation [11]. If we let (p„{x) = exp(^)p„(x), Eq. (l37|) 
becomes 

r J2 

(p„{x) = An(p„{x) , (94) 



where 



D I dWjx) ] _ D (fW{x) 



(95) 



which is exactly the effective potential in Eq.dTst. The eigenvalue spectrum remains the 
same 

We write a double-well potential in the form 



U* =H 



l2 



(96) 



where H is the height of the barrier in units of fcgT and W is the half-width of the barrier. 
We will consider double-wells having high barriers, H » 1, and fixed half width, W. 
The absorbing walls are placed at the two minima, x = ±W. The Schrodinger potential 
corresponding to (l96b is 



V,ix) = D 



6W 



2H UH^ 



8H2 , 4H2 
-X* -I- 



PV6 



PV8 



X 



(97) 



In Fig|9lwe plot U*{x) and Vs{x) for a double-well with sufficiently high barrier. For 
the potential in the Schrodinger picture, Vg = +<x> outside the central interval to ensure 
the wave functions (pn{x) vanish at the ends of the interval, thus satisfying the absorbing 
boundary conditions. 

We use perturbation theory to describe the lowest stationary state, which must exist 
because of the boundary conditions. For sufficiently high barriers, we expect that, the 
lowest eigenstate will be localized at the central minimum, suggesting the use of a pertur- 
bation process based on a simple harmonic oscillator. Using a numerical procedure for 
bound-state solutions of the time-independent Schrodinger equation [73], we confirmed 
this localization for high barriers. We also note that for high barriers, the oscillator's 
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Gaussian wave function nearly vanishes at the boundaries, in approximate satisfaction of 
the proper boundary conditions. 

Our perturbation calculation is therefore based on separating off the dominant har- 
monic component of Vg, noting H » 1, using 



Vs{x) = Vo{x) + V^{x), 



(98) 



where 



Vo{x) 
Vi{x) 



D — + —-x^ ' 



\W2 W4 

6H , 8H2 ^ 4H' 



^1 w^""' W 



(99) 
(100) 



From textbook results for a linear harmonic oscillator, the first eigenvalue and wave 
function are 
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where 



The first order perturbative correction is 



(101) 
(102) 

(103) 



(104) 



which is down by a factor of H from the zero order result. The second order correction is 
down by another factor of H: 



D 



O 



1 



HW2 



(105) 
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The net result for the lowest eigenvalue is thus 



Ai = D 



D 



4H_^ / 1 



(106) 
(107) 



where a = AH/W^ gives the curvature at the top of the barrier. 

Eq. ()107|l shows that the long time behavior is simply linear in the barrier-top curvature 
a, for large values of a and fixed W. This is also the case for the inverted parabolic 
potential, as can be determined from Eq. (l48b , or by performing the same calculation as 
we did for the double-well in this section. In fact, when Dta » 1, Eq. (l48b becomes 



pbit) - 



DaW 



1 / 



1 -erf(WV«72). 



n 




exp 



I exp{-Dat) . 



(108) 



We therefore expect similar linearity with a in the higher barrier limit for any system that 
can be approximated by an inverted parabola and a high order correction. 

We performed numerical checks of the approximation ()107|l . We determined the lowest 
eigenvalue Ai numerically using high moments and Eq. (l63b . We also used path-sampling 
simulations as a consistency check. Fig. [TO] compares the numerical evaluation of Ai 
via Eqs. (l63|t and (|64|) with a path-sampling simulation data for a particular double-well 
potential. 

In Fig. [iTl we compare Eq. (ll07l) with numerical calculations of Ai for double-well 
potentials and inverted parabolae with fixed W and a range of curvatures (10 - 100) and 
plot Ai/Da as a function of curvature «W^. As Eq. (11071) predicts, Ai/Da approaches unity 
for large curvature. 



V. SUMMARY 



We have applied a combination of analytic and numerical techniques to study the 
distribution, pb{t), of the durations of transition-events over a barrier in a one-dimensional 
system undergoing over-damped Langevin dynamics. The typical event duration is much 
shorter than the well-studied first passage time (FPT) [10]; see Figs. [Hand IH The event 
duration scale is the simplest non-trivial measure of the detailed dynamics of an activated 
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process, and we believe it is critical for future quantitative study of dynamics of many 
chemical and biological systems 1^. 

The distribution pi,{t) can be derived from the Fokker-Planck equation with special 
boundary conditions and was subjected to detailed analysis. A number of results are 
obtained, including: (i) the analytic form of the asymptotic short-time behavior {t — > 0), 
which is universal and independent of the potential function; (ii) the first non-universal 
correction to the short-time behavior; (iii) following Gardiner a recursive formulation 
for calculating, exactly, all moments of pi, based solely on the potential function — along 
with approximations for the distribution based on a small number of moments; (iv) a 
high-barrier approximation to the long-time {t — > oo) behavior of pbit)', and (v) a rough 
but simple analytic estimate of the average event duration {i)f„ which generally is sensitive 
to details of the potential. All of the analytic results are confirmed by transition-path- 
sampling simulations. 

A number of interesting questions remain open. Perhaps most centrally, what changes 
in ph can be expected for high-dimensional system with potentially "rough" and complex 
landscapes. A particular interest is in conformational transitions in proteins IzJ]. Fur- 
thermore, it is relevant to consider how non-white noise affect event durations. Finally, 
how can the detailed relationship between ph and first passage-time be quantified, if at 
all? 
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Figure 1: Trajectory for a Brownian particle moving in a double well potential. The left graph is a 
long trajectory with several transition-events. The right graph is the detail of a single transition- 
event cut from the same long trajectory. 




Figure 2: The definition of transition-event duration. In the left graph, two positions are defined 
as the start point a and end point h of the transition. For a transition-event in the right graph, 
transition-event duration is the duration between last time the particle passes a and the first time 
it passes h. 
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Figure 3: Scaled event duration distribution for free diffusion and the linear potential. The left 
graph is for free diffusion with different widths L, where data points are from path-sampling 
simulations and numerical evaluations of Eq.ll38ll. The right graph is for the potential U* = kx 
with different widths L and different slopes k; the points are data from path-sampling simulations 
and the lines are numerical evaluations of Eq. i4lTl . 




Figure 4: Inverted parabolic potential with different boundary conditions. The left graph is the 
inverted parabolic potential U* - -20(1 - x^), with open boundary condition. The right graph is 
the same potential with two absorbing walls at ;c = -1 and x = 1. 
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Figure 5: The event-duration distribution pi, for inverted parabolic potentials with different di- 
mensionless barrier heights ^. The absorbing walls are at -W and W, with W — 1. The data are 
from path-sampling simulation and numerical evaluations (lines) use the approximate formula, 
Eq.lliSl. This approximation improves with increasing barrier height. 
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Figure 6: Reconstruction of pi, for double-well potentials U* - H 
and a Gram-Schmidt procedure. The two absorbing walls are at - W and W, where W = 1. The data 
are from path-sampling simulations and numerical results are based on theory in Section PVD II 
and lIVII)4l using the first five moments. They match well except the long time tail as emphasized 
in the semi-log plot. 
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Figure 7: Calculating the path integral between the start point '+' and the end point 'o' with 
absorbing walls at :c = and x - L. The dashed curve represents an arbitrary potential between 
the two absorbing walls. Eq.llSlTl indicates that one must calculate the path integrals between the 
start point '+' and all "end points", including the real one 'o' and image end points 'A', in the 
periodic potential with open boundary condition. 
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Figure 8: Short time behavior of pj, for double-well potentials U* = H 
sorbing walls at -W and W, where W — 1. Data from path-sampling simulations (solid lines) are 
compared to numerical estimation of pf,(i — > 0) (dashed lines) based on Eqs.iSSl, iSSll and l l9Tl to 
numerically estimate the short time behavior of p^. 




Figure 9: A high barrier double-well potential (left) and its Schrodinger analogue (right) from 
Eqs.i96l and ||97|| , for H - 8 and W — 1. The two absorbing walls are put at the two minima. If the 
barrier height of the potential is not sufficiently large, the minimum of at ^ = will disappear. 
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Figure 10: Comparison of two numerical approaches for computing the long time behavior of 
the event-duration distribution. A path-sampling simulation is compared to pure exponential 
behavior based on the lowest eigenvalue Ai determined from Eqs.i63l and <64t with H = 10 and 
W = 1. Similar results are obtained for other high barrier cases. 
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Figure 11: Long time behavior of pj, for double- well potentials and the inverted parabola. Exact 
numerical results for lowest eigenvalue Ai/Da are plotted as a function of the dimensionless 
curvature aW^ at the barrier peaks. Note that barrier heights are proportional to a for the models 
considered. Two double-well potentials with fixed = 1 and — 3, and one inverted parabola 
system with fixed = 1 are considered. Filled symbols indicate the a values where the barrier 
height H is 10 (in units of hgT). The values of Ai/Da, which dominate the long-time behavior, 
approach 1 for large curvatures, as predicted by Eq. jl07ll . 
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